D-place FARM documentation: Module 3

Module 3: Analysing results produce from Modules 1 and 2

library(png)

Consolidate available files into a single table

# A frequent scenario in our analysis code is that we need to read thousands of 
# output files into a single table. 
# There are a couple flavors of file name that need to be dealt with and the 
# output file needs to include a column with the full file name for reference 
# later. We want to save the final table as it's own file to a folder outside of 
# the one we just read. 

## Here is one flaver where the input is two folders, one with extant 
# simulations and one with extinct simulations. Need to read the files and 
# produce two tables, one for extinct and one for extant. 

## First consolidate the available files into a single table
concatenate <- function(path) {
 available <- list.files(path, full.names = TRUE)
 name <- unlist(strsplit(available[1], split="_"))
 n <- length(available)
 n_name <- length(name) - 1
 load(available[1])
 ncol_files <- length(Sim_statistics[[1]]) + n_name
 files <- matrix(nrow = n,
                 ncol = ncol_files)
 split_names <- function(x){unlist(strsplit(x, split = "_"))}
 name_available <- do.call(rbind, lapply(as.list(available), split_names))
 files[, 1:n_name] <- name_available[, -ncol(name_available)]

 for (i in 1:n) {
   error <- try(load(available[i]), silent = TRUE)
   if (class(error) != "try-error") {
     files[i, (n_name + 1):ncol_files] <- Sim_statistics[[1]]
   }
 }


 files <- files[, -c(1, 2, 3, 5, 7, 8, 13, 18, 23, 28, 33, 35)]
 last_name <- colnames(Sim_statistics[[1]])
 if (is.null(last_name)) {

   files <- files[, -ncol(files)]
 }
 colnames(files) <-  c("Sim_stats_rep", "combo", paste0("P.speciation", 1:4),
                       paste0("P.extinct", 1:4),  paste0("P.diffus", 1:4),
                       paste0("P.TO", 1:4),  paste0("P.Arisal", 1:4), 
                       "timesteps", "NBS", last_name)

 Concatenated_data <- as.data.frame(files)

 begin <- which(colnames(Concatenated_data) ==  "number_of_branches")
 Concatenated_data_stat <- Concatenated_data[, begin:ncol(Concatenated_data)]
 Concatenated_data_stat <- apply(Concatenated_data_stat, 2, as.numeric)
 remove <- which(is.na(rowSums(Concatenated_data_stat)))
 Concatenated_data <- Concatenated_data[-remove, ]

 one <- subset(Concatenated_data, combo =="01")
 two <- subset(Concatenated_data, combo =="02")
 five <- subset(Concatenated_data, combo =="05")
 six <- subset(Concatenated_data, combo =="06")

 crop <- min(sapply(list(one, two, five, six), nrow))

 one <- one[1:crop, ]
 two <- two[1:crop, ]
 five <- five[1:crop, ]
 six <- six[1:crop, ]

 Concatenated_data2 <- rbind(one, two, five, six)
 res <- list(Concatenated_data, Concatenated_data2, crop)
 names(res) <- c("Concatenated_data", "Concatenated_data_crop", "crop")
 return(res)
 
}

path <- "~/Box Sync/Four model compare third run/Module 2"
Concatenated_data0 <- concatenate(path)
print(Concatenated_data0[[3]])
path_save <- "~/Box Sync/Four model compare third run"
Concatenated_data <- Concatenated_data0[[1]]
save(Concatenated_data, file = paste0(path_save,
                                     "/Four_model_compare_results_not_cropped.Rdata"))
Concatenated_data <- Concatenated_data0[[2]]
save(Concatenated_data, file=paste0(path_save,
                                   "/Four_model_compare_results_", 
                                   format(Sys.time(), format="%d_%b_%Y"),
                                   "_crop_to_", Concatenated_data0[[3]],".Rdata"))

### Repeated for extinct
path <- "~/Box Sync/Four model compare third run/Module 2 extinct"
Concatenated_data <- concatenate(path)[[1]]
path_save <- "~/Box Sync/Four model compare third run"

save(Concatenated_data, file=paste0(path_save,
                                   "/Four_model_compare_results_extinct_", 
                                   format(Sys.time(), format="%d_%b_%Y"),".Rdata"))

Summarize overall extinction rates

load("Four_model_compare_results_02_Aug_2017_crop_to_6128.Rdata")
extant <- Concatenated_data
extant

 load("Four_model_compare_results_extinct_02_Aug_2017.Rdata")
extinct <- Concatenated_data
extinct

head(extant)
head(extinct)

for(i in c(3:22)){
    extinct[which(is.nan(as.numeric(as.character(extinct[, i]))) == TRUE), i] <- NA
}

for(i in c(3:22)){
    extant[which(is.nan(as.numeric(as.character(extant[, i]))) == TRUE), i] <- NA
}




xlimit <- c(0,1)
ylimit <- c(0,3000)
maincex <- 0.9

png(file="Global_success_rate_per_parameter.png", width=8.5, height=11, units="in", res=300)

par(mfrow=c(5,4), mar=c(3,3,3,0))


hist(as.numeric(as.character(extinct[,3])), main="speciation of F in F env", col=adjustcolor("firebrick", alpha=0.7), breaks=100, border=NA, xlim= xlimit, ylim= ylimit, cex.main= maincex)
hist(as.numeric(as.character(extant[,3])), main="speciation of F in F env", col=adjustcolor("cornflowerblue", alpha=0.7), breaks=100, border=NA, xlim= xlimit, ylim= ylimit, cex.main= maincex, add=TRUE)


hist(as.numeric(as.character(extinct[,4])), main="speciation of D in F env", col=adjustcolor("firebrick", alpha= 0.7), breaks=100, border=NA, xlim= xlimit, ylim= ylimit, cex.main= maincex)
hist(as.numeric(as.character(extant[,4])), main="speciation of D in F env", col=adjustcolor("cornflowerblue", alpha= 0.7), breaks=100, border=NA, xlim= xlimit, ylim= ylimit, cex.main= maincex, add=TRUE)


hist(as.numeric(as.character(extinct[,5])), main="speciation of F in D env", col=adjustcolor("firebrick", alpha= 0.7), breaks=100, border=NA, xlim= xlimit, ylim= ylimit, cex.main= maincex)
hist(as.numeric(as.character(extant[,5])), main="speciation of F in D env", col=adjustcolor("cornflowerblue", alpha= 0.7), breaks=100, border=NA, xlim= xlimit, ylim= ylimit, cex.main= maincex, add=TRUE)

hist(as.numeric(as.character(extinct[,6])), main="speciation of D in D env", col=adjustcolor("firebrick", alpha= 0.7), breaks=100, border=NA, xlim= xlimit, ylim= ylimit, cex.main= maincex)
hist(as.numeric(as.character(extant[,6])), main="speciation of D in D env", col=adjustcolor("cornflowerblue", alpha= 0.7), breaks=100, border=NA, xlim= xlimit, ylim= ylimit, cex.main= maincex, add=TRUE)

#######

hist(as.numeric(as.character(extinct[, 7])), main="extinction of F in F env", col=adjustcolor("firebrick", alpha= 0.7), breaks=100, border=NA, xlim= xlimit, ylim= ylimit, cex.main= maincex)
hist(as.numeric(as.character(extant[, 7])), main="extinction of F in F env", col=adjustcolor("cornflowerblue", alpha= 0.7), breaks=100, border=NA, xlim= xlimit, ylim= ylimit, cex.main= maincex, add=TRUE)


hist(as.numeric(as.character(extinct[, 8])), main="extinction of D in F env", col=adjustcolor("firebrick", alpha= 0.7), breaks=100, border=NA, xlim= xlimit, ylim= ylimit, cex.main= maincex)
hist(as.numeric(as.character(extant[, 8])), main="extinction of D in F env", col=adjustcolor("cornflowerblue", alpha= 0.7), breaks=100, border=NA, xlim= xlimit, ylim= ylimit, cex.main= maincex, add=TRUE)

hist(as.numeric(as.character(extinct[, 9])), main="extinction of F in D env", col=adjustcolor("firebrick", alpha= 0.7), breaks=100, border=NA, xlim= xlimit, ylim= ylimit, cex.main= maincex)
hist(as.numeric(as.character(extant[, 9])), main="extinction of F in D env", col=adjustcolor("cornflowerblue", alpha= 0.7), breaks=100, border=NA, xlim= xlimit, ylim= ylimit, cex.main= maincex, add=TRUE)

hist(as.numeric(as.character(extinct[, 10])), main="extinction of D in D env", col=adjustcolor("firebrick", alpha= 0.7), breaks=100, border=NA, xlim= xlimit, ylim= ylimit, cex.main= maincex)
hist(as.numeric(as.character(extant[, 10])), main="extinction of D in D env", col=adjustcolor("cornflowerblue", alpha= 0.7), breaks=100, border=NA, xlim= xlimit, ylim= ylimit, cex.main= maincex, add=TRUE)

######

hist(as.numeric(as.character(extinct[, 11])), main="arisal of F in F env", col=adjustcolor("firebrick", alpha= 0.7), breaks=100, border=NA, xlim= xlimit, ylim= ylimit, cex.main= maincex)
hist(as.numeric(as.character(extant[, 11])), main="arisal of F in F env", col=adjustcolor("cornflowerblue", alpha= 0.7), breaks=100, border=NA, xlim= xlimit, ylim= ylimit, cex.main= maincex, add=TRUE)



hist(as.numeric(as.character(extinct[, 12])), main="arisal of D in F env", col=adjustcolor("firebrick", alpha= 0.7), breaks=100, border=NA, xlim= xlimit, ylim= ylimit, cex.main= maincex)
hist(as.numeric(as.character(extant[, 12])), main="arisal of D in F env", col=adjustcolor("cornflowerblue", alpha= 0.7), breaks=100, border=NA, xlim= xlimit, ylim= ylimit, cex.main= maincex, add=TRUE)

hist(as.numeric(as.character(extinct[, 13])), main="arisal of F in D env", col=adjustcolor("firebrick", alpha= 0.7), breaks=100, border=NA, xlim= xlimit, ylim= ylimit, cex.main= maincex)
hist(as.numeric(as.character(extant[, 13])), main="arisal of F in D env", col=adjustcolor("cornflowerblue", alpha= 0.7), breaks=100, border=NA, xlim= xlimit, ylim= ylimit, cex.main= maincex, add=TRUE)

hist(as.numeric(as.character(extinct[, 14])), main="arisal of D in D env", col=adjustcolor("firebrick", alpha= 0.7), breaks=100, border=NA, xlim= xlimit, ylim= ylimit, cex.main= maincex)
hist(as.numeric(as.character(extant[, 14])), main="arisal of D in D env", col=adjustcolor("cornflowerblue", alpha= 0.7), breaks=100, border=NA, xlim= xlimit, ylim= ylimit, cex.main= maincex, add=TRUE)

######

hist(as.numeric(as.character(extinct[, 15])), main="NOPE -- Diffusion: source F, target F", col=adjustcolor("firebrick", alpha= 0.7), breaks=100, border=NA, xlim= xlimit, ylim= c(0,18000), cex.main= maincex)
hist(as.numeric(as.character(extant[, 15])), main="NOPE -- Diffusion: source F, target F", col=adjustcolor("cornflowerblue", alpha= 0.7), breaks=100, border=NA, xlim= xlimit, ylim= c(0,18000), cex.main= maincex, add=TRUE)



hist(as.numeric(as.character(extinct[, 16])), main="Diffusion: source D, target F", col=adjustcolor("firebrick", alpha= 0.7), breaks=100, border=NA, xlim= xlimit, ylim= ylimit, cex.main= maincex)
hist(as.numeric(as.character(extant[, 16])), main="Diffusion: source D, target F", col=adjustcolor("cornflowerblue", alpha= 0.7), breaks=100, border=NA, xlim= xlimit, ylim= ylimit, cex.main= maincex, add=TRUE)

hist(as.numeric(as.character(extinct[, 17])), main="Diffusion: source F, target D", col=adjustcolor("firebrick", alpha= 0.7), breaks=100, border=NA, xlim= xlimit, ylim= ylimit, cex.main= maincex)
hist(as.numeric(as.character(extant[, 17])), main="Diffusion: source F, target D", col=adjustcolor("cornflowerblue", alpha= 0.7), breaks=100, border=NA, xlim= xlimit, ylim= ylimit, cex.main= maincex, add=TRUE)

hist(as.numeric(as.character(extinct[, 18])), main="NOPE -- Diffusion: source D, target D", col=adjustcolor("firebrick", alpha= 0.7), breaks=100, border=NA, xlim= xlimit, ylim= c(0,18000), cex.main= maincex)
hist(as.numeric(as.character(extant[, 18])), main="NOPE -- Diffusion: source D, target D", col=adjustcolor("cornflowerblue", alpha= 0.7), breaks=100, border=NA, xlim= xlimit, ylim= c(0,18000), cex.main= maincex, add=TRUE)

####

hist(as.numeric(as.character(extinct[, 19])), main="Takeover: source F, target F", col=adjustcolor("firebrick", alpha= 0.7), breaks=100, border=NA, xlim= xlimit, ylim= ylimit, cex.main= maincex)
hist(as.numeric(as.character(extant[, 19])), main="Takeover: source F, target F", col=adjustcolor("cornflowerblue", alpha= 0.7), breaks=100, border=NA, xlim= xlimit, ylim= ylimit, cex.main= maincex, add=TRUE)



hist(as.numeric(as.character(extinct[, 20])), main="Takeover: source D, target F", col=adjustcolor("firebrick", alpha= 0.7), breaks=100, border=NA, xlim= xlimit, ylim= ylimit, cex.main= maincex)
hist(as.numeric(as.character(extant[, 20])), main="Takeover: source D, target F", col=adjustcolor("cornflowerblue", alpha= 0.7), breaks=100, border=NA, xlim= xlimit, ylim= ylimit, cex.main= maincex, add=TRUE)

hist(as.numeric(as.character(extinct[, 21])), main="Takeover: source F, target D", col=adjustcolor("firebrick", alpha= 0.7), breaks=100, border=NA, xlim= xlimit, ylim= ylimit, cex.main= maincex)
hist(as.numeric(as.character(extant[, 21])), main="Takeover: source F, target D", col=adjustcolor("cornflowerblue", alpha= 0.7), breaks=100, border=NA, xlim= xlimit, ylim= ylimit, cex.main= maincex, add=TRUE)

hist(as.numeric(as.character(extinct[, 22])), main="Takeover: source D, target D", col=adjustcolor("firebrick", alpha= 0.7), breaks=100, border=NA, xlim= xlimit, ylim= ylimit, cex.main= maincex)
hist(as.numeric(as.character(extant[, 22])), main="Takeover: source D, target D", col=adjustcolor("cornflowerblue", alpha= 0.7), breaks=100, border=NA, xlim= xlimit, ylim= ylimit, cex.main= maincex, add=TRUE)


dev.off()

Run a single random forest on available outputs

# Packages
library(randomForest)

# Radom Forest Function
RF <- function(table_sim, table_real, ntree = 2000, stats_remove = NULL,
              repetitions = 100) {
 # table_sim = table croped 

 begin <- which(colnames(table_sim) ==  "number_of_branches")
 table_sim_data <- table_sim[, begin:ncol(table_sim)]
 table_sim_data <- apply(table_sim_data, 2, as.numeric)
 if (any(is.infinite(table_sim_data))) {
   stop("There is infinite values in the simulation statistics")
 }
 rf_data <- data.frame("Model" = table_sim$combo, table_sim_data)
 rf_data_real <- data.frame(table_real)

 if (!is.null(stats_remove)) {
   rf_data <- rf_data[, -stats_remove]
   rf_data_real <- rf_data_real[, -stats_remove]
 }

 fun <- function(x, y, per = .33) {
   sample(which(y$Model == x), round(table(y$Model)[1]*per))
 }

 results <- matrix(nrow = repetitions, ncol = length(table_real) * 6 + 21)
 for (i in 1:repetitions) {
   sub.test <- unlist(lapply(as.list(paste0(0, c(1,2,5,6))), fun,
                             y = rf_data))
   test2 <- rf_data[sub.test, 2:ncol(rf_data)]
   test1 <- rf_data[sub.test, 1]
   train <- rf_data[-sub.test, ]

   fit <- randomForest(Model ~ ., data = train, xtest = test2, 
                       ytest = test1, importance = TRUE, 
                       ntree = ntree, keep.forest = TRUE,
                       replace = TRUE)
  #dev.new()
  #plot(fit)
   predictions <- predict(fit, rf_data_real, type = "prob")
   var_import <- importance(fit)

   error <- mean(fit$test$confusion[, 5])
   confusion <- as.numeric(fit$test$confusion[, 1:4])
   names(confusion) <- paste0("Confusion_", 
                              rep(colnames(fit$test$confusion[, 1:4]),
                                  each = 4),
                              rownames(fit$test$confusion[, 1:4]))
   predictions_prob <- as.numeric(predictions)
   names(predictions_prob) <- paste0("Prediction_prob", colnames(predictions))
   var_import_vec <- as.numeric(var_import)
   names(var_import_vec) <- paste0("var_imp_", 
                                   rep(colnames(var_import), 
                                       each = nrow(var_import)),
                                   "_",
                                   rownames(var_import))
   vec <- c(error, confusion, predictions_prob, var_import_vec)
   results[i, ] <- vec
 }
 colnames(results) <- names(vec)
 colnames(results)[1] <- "Error_test"
 results <- cbind(1:nrow(results), results)
 colnames(results)[1] <- "Replicate"
 return(results)
}

# Multiple Random Forest
RF_mult <- function(table_sim, table_real_mult, ntree = 2000, 
                   stats_remove = NULL, repetitions = 100) {
 n_r <- nrow(table_real_mult)
 n <- repetitions * n_r
 n_col <- ncol(table_real_mult) * 6 + 22
 results <- matrix(nrow = n, ncol = n_col)
 x <- 0
 for (i in seq(1, n, repetitions)) {
   x <- x + 1
   temp <- RF(table_sim, 
              table_real_mult[x, , drop = FALSE],
              ntree = ntree, 
              repetitions = repetitions, 
              stats_remove = stats_remove)
   results[i:(i + repetitions - 1), ] <- temp
 }
 tree <- rep(1:n_r, each = repetitions)
 results <- cbind(tree, results)
 colnames(results) <- c("Tree", colnames(temp))
 return(results)
}

# Accumulation function
RF_acum <- function(table_sim, table_real_mult, ntree = 2000, 
                   stats_remove = NULL, repetitions = 100,
                   resolution = 100, minimun = 100) {
 sequence <- seq(minimun, nrow(table_sim), resolution)
 n_col <- ncol(table_real_mult) * 6 + 24 
 n_r <- nrow(table_real_mult)
 n_seq <- length(sequence)
 n <- repetitions * n_r * n_seq
 results <- matrix(nrow = n, ncol = n_col)
 lottery1 <- which(table_sim$combo == "01")
 lottery2 <- which(table_sim$combo == "02")
 lottery5 <- which(table_sim$combo == "05")
 lottery6 <- which(table_sim$combo == "06")
 sub1 <- sample(lottery1, sequence[1]/4, replace = FALSE)
 sub2 <- sample(lottery2, sequence[1]/4, replace = FALSE)
 sub5 <- sample(lottery5, sequence[1]/4, replace = FALSE)
 sub6 <- sample(lottery6, sequence[1]/4, replace = FALSE)
 for (i in 1:n_seq) {
   if (i != 1) {
     sub1 <- c(sub1, sample(lottery1[-sub1], sequence[1]/4, replace = FALSE))
     sub2 <- c(sub2, sample(lottery2[-sub2], sequence[1]/4, replace = FALSE))
     sub5 <- c(sub5, sample(lottery5[-sub5], sequence[1]/4, replace = FALSE))
     sub6 <- c(sub6, sample(lottery6[-sub6], sequence[1]/4, replace = FALSE))
   }
   temp <- RF_mult(table_sim[c(sub1, sub2, sub5, sub6), ], table_real_mult,
                   ntree = ntree, repetitions = repetitions,
                   stats_remove = stats_remove)
   n_temp <- nrow(temp)
   results[1:n_temp + (n_temp * (i - 1)), ] <- cbind(rep(sequence[i],
                                                         n_temp), temp)
 }
 colnames(results) <- c("Subsample_size", colnames(temp))
 return(results)
}

# TEST FINAL (RUN ONLY THIS)
# Load data
setwd("~/Desktop/RF after bruno")
load("Four model compare third runFour_model_compare_results_28_Jul_2017_crop_to_5767.Rdata")
load("real.analysis_mult.RData")

table_sim = Concatenated_data
table_real = real.analysis.mult[[1]][[1]]
table_real_mult = do.call(rbind, lapply(real.analysis.mult, 
                                       function(x){x[[1]]}))

#start time
time_start <- Sys.time()

# Change the parameters as you wish
res_acum <- RF_acum(table_sim, table_real_mult[100:110,,drop=FALSE], ntree = 2000, 
                   stats_remove = NULL, repetitions = 1,
                   resolution = 5000, minimun = 1000)
                   
#stop time
time_stop <- Sys.time()
difftime(time_stop, time_start)                   

save(res_acum, file="~/Desktop/10_more_tree_out.Rdata")

# Change argument FUN to sd to obtain the standard deviation
plot(aggregate(x = res_acum[, 4], by = list((res_acum[, 1])), FUN = mean), type = "n", ylim=c(0,1))

lines(aggregate(x = res_acum[, 21], by = list((res_acum[, 1])), FUN = mean), type = "l", col="orange")
lines(aggregate(x = res_acum[, 22], by = list((res_acum[, 1])), FUN = mean), type = "l", col="blue")
lines(aggregate(x = res_acum[, 23], by = list((res_acum[, 1])), FUN = mean), type = "l", col="pink")
lines(aggregate(x = res_acum[, 24], by = list((res_acum[, 1])), FUN = mean), type = "l", col="darkgreen")

Visualize outputs from Random Forest analysis

# load outputs from RF runs
load("/Users/Ty/Desktop/small_RF_out.Rdata")
a <- res_acum

load("/Users/Ty/Desktop/10_more_tree_out.Rdata")
b <- res_acum

load("/Users/Ty/Desktop/10_tree_out.Rdata")
c <- res_acum

#object called res_acum
res_acum <- rbind(b,c)

png("overall_error_per_sample_size.png", width = 11, height = 8.5, res = 300, units = "in")
plot(aggregate(x = res_acum[, 4], by = list((res_acum[, 1])), FUN = mean), type="l", ylim=c(0,0.5), ylab="overall confusion error", xlab="sample size")
dev.off()

null device 
          1 

png("predictions.png", width = 11, height = 8.5, res = 300, units = "in")
par(mar=c(5,6,2,1))
# Change argument FUN to sd to obtain the standard deviation
plot(aggregate(x = res_acum[, 4], by = list((res_acum[, 1])), FUN = mean), type = "n", ylim=c(0,1), ylab="probability that our known cultural phylogeny came \n from each type of simulated mechanism", xlab="number of simulation replicates", cex.axis=1.5, cex.lab=1.5)
basic_mean <- aggregate(x = res_acum[, 21], by = list((res_acum[, 1])), FUN = mean)
diffusion_mean <- aggregate(x = res_acum[, 22], by = list((res_acum[, 1])), FUN = mean)
TO_mean <- aggregate(x = res_acum[, 23], by = list((res_acum[, 1])), FUN = mean)
both_mean <- aggregate(x = res_acum[, 24], by = list((res_acum[, 1])), FUN = mean)
                   
basic_sd <- aggregate(x = res_acum[, 21], by = list((res_acum[, 1])), FUN = sd)
diffusion_sd <- aggregate(x = res_acum[, 22], by = list((res_acum[, 1])), FUN = sd)
TO_sd <- aggregate(x = res_acum[, 23], by = list((res_acum[, 1])), FUN = sd)
both_sd <- aggregate(x = res_acum[, 24], by = list((res_acum[, 1])), FUN = sd)
polygon(x=c(basic_mean[,1], rev(basic_mean[,1])), y=c(basic_mean[,2] + basic_sd[,2], rev(basic_mean[,2] - basic_sd[,2])), col=adjustcolor("orange", alpha=0.5), border=NA)
polygon(x=c(diffusion_mean[,1], rev(diffusion_mean[,1])), y=c(diffusion_mean[,2] + diffusion_sd[,2], rev(diffusion_mean[,2] - diffusion_sd[,2])), col=adjustcolor("cornflowerblue", alpha=0.5), border=NA)
polygon(x=c(TO_mean[,1], rev(TO_mean[,1])), y=c(TO_mean[,2] + TO_sd[,2], rev(TO_mean[,2] - TO_sd[,2])), col=adjustcolor("pink", alpha=0.5), border=NA)
polygon(x=c(both_mean[,1], rev(both_mean[,1])), y=c(both_mean[,2] + both_sd[,2], rev(both_mean[,2] - both_sd[,2])), col=adjustcolor("darkgreen", alpha=0.5), border=NA)
lines(basic_mean,  col="orange")
lines(diffusion_mean,  col="cornflowerblue")
lines(TO_mean,  col="pink")
lines(both_mean,  col="darkgreen")
labs <- c("Basic", "+Diffusion", "+Takeover", "+Diffusion +Takeover")
legend("topright", legend=labs, col=c("orange", "cornflowerblue", "pink", "darkgreen"), lty=1, cex=1.5, lwd=10)
dev.off()

null device 
          1 

png("last_step_predictions.png", width = 11, height = 8.5, res = 300, units = "in")
par(mar=c(5,6,2,1))
this <- length(basic_mean[,1])

plot(x=c(0,5), y=c(0,1), type="n",  xaxt="n", ylab="probability that our known cultural phylogeny came \n from each type of simulated mechanism", xlab="", xlim=c(0.5,4.5))
polygon(x=c(0.75,0.75,1.25,1.25), y=c(basic_mean[this,2] - basic_sd[this,2], basic_mean[this,2] + basic_sd[this,2], basic_mean[this,2] + basic_sd[this,2], basic_mean[this,2] - basic_sd[this,2]), col = adjustcolor("orange", alpha=0.5), border=NA )
polygon(x=c(1.75,1.75,2.25,2.25), y=c(diffusion_mean[this,2] - diffusion_sd[this,2], diffusion_mean[this,2] + diffusion_sd[this,2], diffusion_mean[this,2] + diffusion_sd[this,2], diffusion_mean[this,2] - diffusion_sd[this,2]), col = adjustcolor("cornflowerblue", alpha=0.5), border=NA)
polygon(x=c(2.75,2.75,3.25,3.25), y=c(TO_mean[this,2] - TO_sd[this,2], TO_mean[this,2] + TO_sd[this,2], TO_mean[this,2] + TO_sd[this,2], TO_mean[this,2] - TO_sd[this,2]), col = adjustcolor("pink", alpha=0.5), border=NA)
polygon(x=c(3.75,3.75,4.25,4.25), y=c(both_mean[this,2] - both_sd[this,2], both_mean[this,2] + both_sd[this,2], both_mean[this,2] + both_sd[this,2], both_mean[this,2] - both_sd[this,2]), col = adjustcolor("darkgreen", alpha=0.5), border=NA)


lines(x=c(0.5, 1.5), c(basic_mean[this,2], basic_mean[this,2]),  col="orange")
lines(x=c(1.5, 2.5), c(diffusion_mean[this,2], diffusion_mean[this,2]),  col="cornflowerblue")
lines(x=c(2.5, 3.5), c(TO_mean[this,2], TO_mean[this,2]),  col="pink")
lines(x=c(3.5, 4.5), c(both_mean[this,2], both_mean[this,2]),  col="darkgreen")

legend("topright", legend=labs, col=c("orange", "cornflowerblue", "pink", "darkgreen"), lty=1, cex=1.5, lwd=10)

dev.off()

#colnames(res_acum)
long <- length(Confusion_0101_mean[,1])
Confusion_sum_01 <- aggregate(x = res_acum[, 5:8], by = list((res_acum[, 1])), FUN = mean)
Confusion_sum_02 <- aggregate(x = res_acum[, 9:12], by = list((res_acum[, 1])), FUN = mean)
Confusion_sum_03 <- aggregate(x = res_acum[, 13:16], by = list((res_acum[, 1])), FUN = mean)
Confusion_sum_04 <- aggregate(x = res_acum[, 17:20], by = list((res_acum[, 1])), FUN = mean)
percent_one <- (Confusion_sum_01[long,2:5]/sum(Confusion_sum_01[long,2:5])) * 100
percent_two <- (Confusion_sum_02[long,2:5]/sum(Confusion_sum_02[long,2:5])) * 100
percent_three <- (Confusion_sum_03[long,2:5]/sum(Confusion_sum_03[long,2:5])) * 100
percent_four <- (Confusion_sum_04[long,2:5]/sum(Confusion_sum_04[long,2:5])) * 100
confusion <- matrix(c(percent_one, percent_two, percent_three, percent_four), 4,4)
labs <- c("Basic", "+Diffusion", "+Takeover", "+Diffusion",  "+Takeover")
colors1 <- colorRampPalette(colors = c("grey95", "grey40"))
png("confusion_matrix.png", width = 8.5, height = 8.5, res = 600, units = "in")
par(mar=c(8,8,1,1))
plot(0,0,xlim=c(-0.2,1.4), ylim=c(-0.2,1.4), xaxt="n", xlab="", yaxt="n", ylab="" , bty="n")
#image(prop, col = colors1(20), axes=FALSE)
axis(1, at=c(0, .4, .8, 1.2, 1.3), labels=labs, tick = FALSE, line = FALSE, cex.axis = 1, pos = -.19, las=2)
axis(2, at=rev(c(-0.05, 0.05, .4, .8, 1.2)), labels=labs, tick = FALSE, line = FALSE, cex.axis = 1, las=2)
mtext("percent of time that RF identifies input model as each model type", side = 1, padj = 10, cex = 1)
mtext("known model type given to random forest", side = 2, padj = -10, cex = 1)
for(i in 1:4) {
  for(j in 4:1) {
    xs <- c(0, .4, .8, 1.2)[i]
    ys <- rev(c(0, .4, .8, 1.2))[j]
    polygon(x=c(xs-0.2, xs-0.2, xs+0.2, xs+0.2), y=c(ys-0.2, ys+0.2, ys+0.2, ys-0.2), col=colors1(100)[round(as.numeric(confusion[i, j]), 1)+1])
    if(i == j){text(x = xs, y = ys, paste0(round(as.numeric(confusion[i, j]), 2), "%"), cex = 2.2, col="limegreen")}else{(text(x = xs, y = ys, paste0(round(as.numeric(confusion[i, j]), 2), "%"), cex = 2.2, col="black"))}
  }
}

importance(fit)

# Variables importance

imp <- importance(fit)
imp <- apply(imp, 2, function(x) (x - min(x))/(max(x) - min(x)))
imp <- imp[sort(imp[, 5], index.return = TRUE, decreasing = TRUE)$ix, ]


names <- rownames(imp)
names[names == "spatial.tests.fora"] <- "Space F"
names[names == "spatial.tests.dom"] <- "Space D"
names[names == "sprate"] <- "Sp(ratio)"
names[names == "transition_from_trait_1_to_2"] <- "TR(1-2)"
names[names == "transition_from_trait_2_to_1"] <- "TR(2-1)"
names[names == "Phylogenetic_signal"] <- "PhySig(D)"
names[names == "Evolutionary_distinctiveness_sum"] <- "EDsum"
names[names == "Pylo_diversity_is_sum_of_BL"] <- "PDsum"
names[names == "transition_rate_ratio_1to2_over_2to1"] <- "TR(ratio)"
names[names == "gamma"] <- "Gamma"
names[names == "mean_Phylogenetic_isolation"] <- "MPI"
names[names == "extrate"] <- "Ext(ratio)"
names[names == "average_phylogenetic_diversity_is_mean_of_BL"] <- "PDmean"
names[names == "extinction_per_speciation"] <- "DR"
names[names == "variance_Phylogenetic_isolation"] <- "VPI"
names[names == "F_quadratic_entropy_is_sum_of_PD"] <- "F"
names[names == "Mean_pairwise_distance"] <- "MPD"
names[names == "variance_Pylo_diversity_is_variance_of_BL"] <- "PDvar"
names[names == "variance_pairwise_distance"] <- "VPD"


png("var_import_all.png", width = 25, height = 25, unit="in", res=300)
par(mar = c(10, 18, 1, 1))
plot(x = rev(imp[, 5]), y = 1:nrow(imp), type = "l", yaxt = "n", 
     ylab = "", xlab = "Variable Importance",
     xlim = c(0, 1), lwd = 2, cex.lab = 4)
for (i in 1:nrow(imp)) {
  abline(h = i, lty = 3, col = "gray80")
}
abline(v = seq(0, 1, 1/19), lty = 3, col = "gray80")

lines(x = rev(imp[, 4]), y = 1:nrow(imp), col = "darkgreen", lwd = 2)
lines(x = rev(imp[, 3]), y = 1:nrow(imp), col = "red", lwd = 2)
lines(x = rev(imp[, 2]), y = 1:nrow(imp), col = "blue", lwd = 2)
lines(x = rev(imp[, 1]), y = 1:nrow(imp), col = "darkorange1", lwd = 2)
lines(x = rev(imp[, 5]), y = 1:nrow(imp), lwd = 3)

points(x = rev(imp[, 4]), y = 1:nrow(imp), col = "darkgreen", cex = 2)
points(x = rev(imp[, 3]), y = 1:nrow(imp), col = "red", cex = 2)
points(x = rev(imp[, 2]), y = 1:nrow(imp), col = "blue", cex = 2)
points(x = rev(imp[, 1]), y = 1:nrow(imp), col = "darkorange1", cex = 2)
points(x = rev(imp[, 5]), y = 1:nrow(imp), pch = 20, cex = 3)


text(y = 1:nrow(imp), x = par("usr")[1] - .17, labels = rev(names),
     srt = 0, pos = 4, xpd = T, cex = 4)
dev.off()

par(mfrow=c(2,3))

# Box plots
boxplot(spatial.tests.fora ~ Model, data = data.analysis.comp3)
abline(h = a$spatial.tests.fora, col = "red", lty = 2)

boxplot(spatial.tests.dom ~ Model, data = data.analysis.comp3)
abline(h = a$spatial.tests.fora, col = "red", lty = 2)

boxplot(log(sprate) ~ Model, data = data.analysis.comp3, ylim = c(-10, 10))
abline(h = log(a$sprate), col = "red", lty = 2)

boxplot(log(extrate) ~ Model, data = data.analysis.comp3, ylim = c(-10, 10))
abline(h = log(a$extrate), col = "red", lty = 2)

boxplot(log(transition_rate_ratio_1to2_over_2to1) ~ Model, data = data.analysis.comp3)
abline(h = log(a$sprate), col = "red", lty = 2)

boxplot(Phylogenetic_signal ~ Model, data = data.analysis.comp3, ylim = c(0, 1))
abline(h = a$Phylogenetic_signal, col = "red", lty = 2)

#build a data tracking table to track parameter changes through time

str(fit)

y